rm(list = ls())
load("dataBJPOLS.RData")

inst.1 <- "judgeiv_hd"
endo.1 <- "pti"
outc.1 <- "vote_post"

time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow)"
demo.controls <- "age + I(age^2) +  as.factor(race4) + as.factor(race4) + female + vote_pre + as.factor(noteli) + regis_before"
case.controls <- "as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)"

form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

m1a1 <- ivreg(form.1, data = last.cases)
m1a2 <- ivreg(form.2, data = last.cases)
m1a3 <- ivreg(form.3, data = last.cases)

last.cases$comp_w <- NA

for(m in 0:1) {
  for(j in 0:1) {
    form.3.cc <- formula(paste(endo.1, "~", inst.1, "+", time.controls))
    out <- lm(form.3.cc, data = last.cases[any_prior_case == j & severity == m, ])  
    quants.iv <- quantile(last.cases$judgeiv_hd[last.cases$any_prior_case == j  & last.cases$severity == m], c(0.01, 0.99))
    frac_complier <- out$coefficients["judgeiv_hd"] * quants.iv[2] - out$coefficients["judgeiv_hd"] * quants.iv[1]
    mean_frac <- sum(last.cases$any_prior_case == j  & last.cases$severity == m)/nrow(last.cases)
    last.cases$comp_w[last.cases$any_prior_case == j  & last.cases$severity == m] <- frac_complier/mean_frac
  }
}

## -------------------------------
## Fraction of compliers:
## we define compliers as defendants whose pre-trial detention
## decision would have been different had their case been assigned to the most strict instead of the
## most lenient judge
## -------------------------------

out.1 <- lm(form.3.cc, data = last.cases)  
quants.iv <- quantile(last.cases$judgeiv_hd, c(0.01, 0.99))
coef(out.1)[2] * quants.iv[2] - coef(out.1)[2] * quants.iv[1]


## -------------------------------
## Always Takers:
## Always takers are defendants who would always be detained before trial regardless of the bail
## judge assigned to their case.
## -------------------------------

out.1 <- lm(form.3.cc, data = last.cases)  
quants.iv <- quantile(last.cases$judgeiv_hd, c(0.01, 0.99))
coef(out.1)[1] + coef(out.1)[2] * quants.iv[1]

## -------------------------------
## Never Takers:
## Are defendants who would never be detained before trial
## -------------------------------

out.1 <- lm(form.3.cc, data = last.cases)  
quants.iv <- quantile(last.cases$judgeiv_hd, c(0.01, 0.99))
1 - coef(out.1)[1] - coef(out.1)[2] * quants.iv[2]

w.mean <- function(x, y) { sum(x * y,na.rm=T)/sum(y,na.rm=T) }
se.w.mean <- function(x, y) { sqrt(var(x * y,na.rm=T))/sqrt(sum(y,na.rm=T)) }

## Turnout:
c1.1 <- mean(last.cases$vote_post); c1.2 <- sqrt(var(last.cases$vote_post))/sqrt(nrow(last.cases))
c1.3 <- w.mean(last.cases$vote_post, last.cases$comp_w); c1.4 <- se.w.mean(last.cases$vote_post, last.cases$comp_w)

## Age
c2.1 <- mean(last.cases$age); c2.2 <- sqrt(var(last.cases$age))/sqrt(nrow(last.cases))
c2.3 <- w.mean(last.cases$age, last.cases$comp_w); c2.4 <- se.w.mean(last.cases$age, last.cases$comp_w)

## Gender
c3.1 <- mean(last.cases$female); c3.2 <- sqrt(var(last.cases$female))/sqrt(nrow(last.cases))
c3.3 <- w.mean(last.cases$female, last.cases$comp_w); c3.4 <- se.w.mean(last.cases$female, last.cases$comp_w)

## Black
c4.1 <- mean(last.cases$race4 == "B"); c4.2 <- sqrt(var(last.cases$race4 == "B"))/sqrt(nrow(last.cases))
c4.3 <- w.mean(last.cases$race4 == "B", last.cases$comp_w); c4.4 <- se.w.mean(last.cases$race4 == "B", last.cases$comp_w)

## White
c5.1 <- mean(last.cases$race4 == "W"); c5.2 <- sqrt(var(last.cases$race4 == "W"))/sqrt(nrow(last.cases))
c5.3 <- w.mean(last.cases$race4 == "W", last.cases$comp_w); c5.4 <- se.w.mean(last.cases$race4 == "W", last.cases$comp_w)

## Hispanics
c6.1 <- mean(last.cases$race4 == "H"); c6.2 <- sqrt(var(last.cases$race4 == "H"))/sqrt(nrow(last.cases))
c6.3 <- w.mean(last.cases$race4 == "H", last.cases$comp_w); c6.4 <- se.w.mean(last.cases$race4 == "H", last.cases$comp_w)

## any_drug
c7.1 <- mean(last.cases$any_drug); c7.2 <- sqrt(var(last.cases$any_drug))/sqrt(nrow(last.cases))
c7.3 <- w.mean(last.cases$any_drug, last.cases$comp_w); c7.4 <- se.w.mean(last.cases$any_drug, last.cases$comp_w)

## fire_arms_2
c9.1 <- mean(last.cases$any_weapon); c9.2 <- sqrt(var(last.cases$any_weapon))/sqrt(nrow(last.cases))
c9.3 <- w.mean(last.cases$any_weapon, last.cases$comp_w); c9.4 <- se.w.mean(last.cases$any_weapon, last.cases$comp_w)

## fire_arms_2
c10.1 <- mean(last.cases$any_prop); c10.2 <- sqrt(var(last.cases$any_prop))/sqrt(nrow(last.cases))
c10.3 <- w.mean(last.cases$any_prop, last.cases$comp_w); c10.4 <- se.w.mean(last.cases$any_prop, last.cases$comp_w)

## prior_2
c11.1 <- mean(last.cases$any_prior_case); c11.2 <- sqrt(var(last.cases$any_prior_case))/sqrt(nrow(last.cases))
c11.3 <- w.mean(last.cases$any_prior_case, last.cases$comp_w); c11.4 <- se.w.mean(last.cases$any_prior_case, last.cases$comp_w)

## voted 2008
c12.1 <- mean(last.cases$vote_pre); c12.2 <- sqrt(var(last.cases$vote_pre))/sqrt(nrow(last.cases))
c12.3 <- w.mean(last.cases$vote_pre, last.cases$comp_w); c12.4 <- se.w.mean(last.cases$vote_pre, last.cases$comp_w)

## not eligible 2008
c13.1 <- mean(last.cases$noteli); c13.2 <- sqrt(var(last.cases$noteli))/sqrt(nrow(last.cases))
c13.3 <- w.mean(last.cases$noteli, last.cases$comp_w); c13.4 <- se.w.mean(last.cases$noteli, last.cases$comp_w)

## Income Below
c14.1 <- mean(last.cases$median_income_group == "Below"); c14.2 <- sqrt(var(last.cases$median_income_group == "Below"))/sqrt(nrow(last.cases))
c14.3 <- w.mean(last.cases$median_income_group == "Below", last.cases$comp_w); c14.4 <- se.w.mean(last.cases$median_income_group == "Below", last.cases$comp_w)

## Income Above
c15.1 <- mean(last.cases$median_income_group == "Above"); c15.2 <- sqrt(var(last.cases$median_income_group == "Above"))/sqrt(nrow(last.cases))
c15.3 <- w.mean(last.cases$median_income_group == "Above", last.cases$comp_w); c15.4 <- se.w.mean(last.cases$median_income_group == "Above", last.cases$comp_w)

## Income Not Avail
c16.1 <- mean(last.cases$median_income_group == "Unavailable"); c16.2 <- sqrt(var(last.cases$median_income_group == "Unavailable"))/sqrt(nrow(last.cases))
c16.3 <- w.mean(last.cases$median_income_group == "Unavailable", last.cases$comp_w); c16.4 <- se.w.mean(last.cases$median_income_group == "Unavailable", last.cases$comp_w)

tab <- rbind(cbind(c1.1, c1.2, c1.3, c1.4),
             cbind(c2.1, c2.2, c2.3, c2.4),
             cbind(c3.1, c3.2, c3.3, c3.4),
             cbind(c4.1, c4.2, c4.3, c4.4),
             cbind(c5.1, c5.2, c5.3, c5.4),
             cbind(c6.1, c6.2, c6.3, c6.4),
             cbind(c7.1, c7.2, c7.3, c7.4),
             cbind(c9.1, c9.2, c9.3, c9.4),
             cbind(c10.1, c10.2, c10.3, c10.4),
             cbind(c11.1, c11.2, c11.3, c11.4),
             cbind(c12.1, c12.2, c12.3, c12.4),
             cbind(c13.1, c13.2, c13.3, c13.4),
             cbind(c14.1, c14.2, c14.3, c14.4),
             cbind(c15.1, c15.2, c15.3, c15.4),
             cbind(c16.1, c16.2, c16.3, c16.4)
)

colnames(tab) <- c("Mean", "s.e.", "Weighted Mean", "s.e")
rownames(tab) <- c("Turnout", "Age", "Female", "Black", "White", "Hispanic",
                   "Drug-related", "Weapons", "Robbery", 
                   "Prior Offense", "Previous Turnout", "Non-Eligibility",
                   "Income Below", "Income Above", "Income Unk")

tab.o <- data.table(t(tab[, c(1, 3)]))

run <- T
if(run == T) {
  library(doParallel)
  nt <- detectCores()
  nc <- makeCluster(nt)
  registerDoParallel(nc)
  
  results <- foreach(i = 1:500) %dopar% {
    library('AER')
    library('data.table')
    
    out <- last.cases[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
    last.cases.b <- last.cases[last.cases$this_id_fa %in% out$V1, ]
    
    inst.1 <- "judgeiv_hd"
    endo.1 <- "pti"
    outc.1 <- "vote_post"
    
    time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow)"
    demo.controls <- "age + I(age^2) +  as.factor(race4) + as.factor(race4) + female + vote_pre + as.factor(noteli) + regis_before"
    case.controls <- "as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)"
    
    form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
    form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
    form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))
    
    m1a1 <- ivreg(form.1, data = last.cases.b)
    m1a2 <- ivreg(form.2, data = last.cases.b)
    m1a3 <- ivreg(form.3, data = last.cases.b)
    
    last.cases.b$comp_w <- NA
    for(m in 0:1) {
      for(j in 0:1) {
        form.3.cc <- formula(paste(endo.1, "~", inst.1, "+", time.controls))
        out <- lm(form.3.cc, data = last.cases.b[any_prior_case == j & severity == m, ])  
        quants.iv <- quantile(last.cases.b$judgeiv_hd[last.cases.b$any_prior_case == j  & last.cases.b$severity == m], c(0.01, 0.99))
        frac_complier <- out$coefficients["judgeiv_hd"] * quants.iv[2] - out$coefficients["judgeiv_hd"] * quants.iv[1]
        mean_frac <- sum(last.cases.b$any_prior_case == j  & last.cases.b$severity == m)/nrow(last.cases.b)
        last.cases.b$comp_w[last.cases.b$any_prior_case == j  & last.cases.b$severity == m] <- frac_complier/mean_frac
      }
    }
    
    
    ## Fraction of compliers:
    ## we define compliers as defendants whose pre-trial detention
    ## decision would have been different had their case been assigned to the most strict instead of the
    ## most lenient judge
    ## -------------------------------
    
    out.1 <- lm(form.3.cc, data = last.cases.b)  
    quants.iv <- quantile(last.cases.b$judgeiv_hd, c(0.01, 0.99))
    coef(out.1)[2] * quants.iv[2] - coef(out.1)[2] * quants.iv[1]
    
    
    ## -------------------------------
    ## Always Takers:
    ## Always takers are defendants who would always be detained before trial regardless of the bail
    ## judge assigned to their case.
    ## -------------------------------
    
    out.1 <- lm(form.3.cc, data = last.cases.b)  
    quants.iv <- quantile(last.cases.b$judgeiv_hd, c(0.01, 0.99))
    coef(out.1)[1] + coef(out.1)[2] * quants.iv[1]
    
    ## -------------------------------
    ## Never Takers:
    ## Are defendants who would never be detained before trial
    ## -------------------------------
    
    out.1 <- lm(form.3.cc, data = last.cases.b)  
    quants.iv <- quantile(last.cases.b$judgeiv_hd, c(0.01, 0.99))
    1 - coef(out.1)[1] - coef(out.1)[2] * quants.iv[2]
    
    w.mean <- function(x, y) { sum(x * y,na.rm=T)/sum(y,na.rm=T) }
    se.w.mean <- function(x, y) { sqrt(var(x * y,na.rm=T))/sqrt(sum(y,na.rm=T)) }
    
    ## Turnout:
    c1.1 <- mean(last.cases.b$vote_post); c1.2 <- sqrt(var(last.cases.b$vote_post))/sqrt(nrow(last.cases.b))
    c1.3 <- w.mean(last.cases.b$vote_post, last.cases.b$comp_w); c1.4 <- se.w.mean(last.cases.b$vote_post, last.cases.b$comp_w)
    
    ## Age
    c2.1 <- mean(last.cases.b$age); c2.2 <- sqrt(var(last.cases.b$age))/sqrt(nrow(last.cases.b))
    c2.3 <- w.mean(last.cases.b$age, last.cases.b$comp_w); c2.4 <- se.w.mean(last.cases.b$age, last.cases.b$comp_w)
    
    ## Gender
    c3.1 <- mean(last.cases.b$female); c3.2 <- sqrt(var(last.cases.b$female))/sqrt(nrow(last.cases.b))
    c3.3 <- w.mean(last.cases.b$female, last.cases.b$comp_w); c3.4 <- se.w.mean(last.cases.b$female, last.cases.b$comp_w)
    
    ## Black
    c4.1 <- mean(last.cases.b$race4 == "B"); c4.2 <- sqrt(var(last.cases.b$race4 == "B"))/sqrt(nrow(last.cases.b))
    c4.3 <- w.mean(last.cases.b$race4 == "B", last.cases.b$comp_w); c4.4 <- se.w.mean(last.cases.b$race4 == "B", last.cases.b$comp_w)
    
    ## White
    c5.1 <- mean(last.cases.b$race4 == "W"); c5.2 <- sqrt(var(last.cases.b$race4 == "W"))/sqrt(nrow(last.cases.b))
    c5.3 <- w.mean(last.cases.b$race4 == "W", last.cases.b$comp_w); c5.4 <- se.w.mean(last.cases.b$race4 == "W", last.cases.b$comp_w)
    
    ## Hispanics
    c6.1 <- mean(last.cases.b$race4 == "H"); c6.2 <- sqrt(var(last.cases.b$race4 == "H"))/sqrt(nrow(last.cases.b))
    c6.3 <- w.mean(last.cases.b$race4 == "H", last.cases.b$comp_w); c6.4 <- se.w.mean(last.cases.b$race4 == "H", last.cases.b$comp_w)
    
    ## any_drug
    c7.1 <- mean(last.cases.b$any_drug); c7.2 <- sqrt(var(last.cases.b$any_drug))/sqrt(nrow(last.cases.b))
    c7.3 <- w.mean(last.cases.b$any_drug, last.cases.b$comp_w); c7.4 <- se.w.mean(last.cases.b$any_drug, last.cases.b$comp_w)
    
    ## fire_arms_2
    c9.1 <- mean(last.cases.b$any_weapon); c9.2 <- sqrt(var(last.cases.b$any_weapon))/sqrt(nrow(last.cases.b))
    c9.3 <- w.mean(last.cases.b$any_weapon, last.cases.b$comp_w); c9.4 <- se.w.mean(last.cases.b$any_weapon, last.cases.b$comp_w)
    
    ## fire_arms_2
    c10.1 <- mean(last.cases.b$any_prop); c10.2 <- sqrt(var(last.cases.b$any_prop))/sqrt(nrow(last.cases.b))
    c10.3 <- w.mean(last.cases.b$any_prop, last.cases.b$comp_w); c10.4 <- se.w.mean(last.cases.b$any_prop, last.cases.b$comp_w)
    
    ## prior_2
    c11.1 <- mean(last.cases.b$any_prior_case); c11.2 <- sqrt(var(last.cases.b$any_prior_case))/sqrt(nrow(last.cases.b))
    c11.3 <- w.mean(last.cases.b$any_prior_case, last.cases.b$comp_w); c11.4 <- se.w.mean(last.cases.b$any_prior_case, last.cases.b$comp_w)
    
    ## voted 2008
    c12.1 <- mean(last.cases.b$vote_pre); c12.2 <- sqrt(var(last.cases.b$vote_pre))/sqrt(nrow(last.cases.b))
    c12.3 <- w.mean(last.cases.b$vote_pre, last.cases.b$comp_w); c12.4 <- se.w.mean(last.cases.b$vote_pre, last.cases.b$comp_w)
    
    ## not eligible 2008
    c13.1 <- mean(last.cases.b$noteli); c13.2 <- sqrt(var(last.cases.b$noteli))/sqrt(nrow(last.cases.b))
    c13.3 <- w.mean(last.cases.b$noteli, last.cases.b$comp_w); c13.4 <- se.w.mean(last.cases.b$noteli, last.cases.b$comp_w)
    
    ## Income Below
    c14.1 <- mean(last.cases.b$median_income_group == "Below"); c14.2 <- sqrt(var(last.cases.b$median_income_group == "Below"))/sqrt(nrow(last.cases.b))
    c14.3 <- w.mean(last.cases.b$median_income_group == "Below", last.cases.b$comp_w); c14.4 <- se.w.mean(last.cases.b$median_income_group == "Below", last.cases.b$comp_w)
    
    ## Income Above
    c15.1 <- mean(last.cases.b$median_income_group == "Above"); c15.2 <- sqrt(var(last.cases.b$median_income_group == "Above"))/sqrt(nrow(last.cases.b))
    c15.3 <- w.mean(last.cases.b$median_income_group == "Above", last.cases.b$comp_w); c15.4 <- se.w.mean(last.cases.b$median_income_group == "Above", last.cases.b$comp_w)
    
    ## Income Not Avail
    c16.1 <- mean(last.cases.b$median_income_group == "Unavailable"); c16.2 <- sqrt(var(last.cases.b$median_income_group == "Unavailable"))/sqrt(nrow(last.cases.b))
    c16.3 <- w.mean(last.cases.b$median_income_group == "Unavailable", last.cases.b$comp_w); c16.4 <- se.w.mean(last.cases.b$median_income_group == "Unavailable", last.cases.b$comp_w)
    
    tab <- rbind(cbind(c1.1, c1.2, c1.3, c1.4),
                 cbind(c2.1, c2.2, c2.3, c2.4),
                 cbind(c3.1, c3.2, c3.3, c3.4),
                 cbind(c4.1, c4.2, c4.3, c4.4),
                 cbind(c5.1, c5.2, c5.3, c5.4),
                 cbind(c6.1, c6.2, c6.3, c6.4),
                 cbind(c7.1, c7.2, c7.3, c7.4),
                 cbind(c9.1, c9.2, c9.3, c9.4),
                 cbind(c10.1, c10.2, c10.3, c10.4),
                 cbind(c11.1, c11.2, c11.3, c11.4),
                 cbind(c12.1, c12.2, c12.3, c12.4),
                 cbind(c13.1, c13.2, c13.3, c13.4),
                 cbind(c14.1, c14.2, c14.3, c14.4),
                 cbind(c15.1, c15.2, c15.3, c15.4),
                 cbind(c16.1, c16.2, c16.3, c16.4)
    )
    
    colnames(tab) <- c("Mean", "s.e.", "Weighted Mean", "s.e")
    rownames(tab) <- c("Turnout", "Age", "Female", "Black", "White", "Hispanic",
                       "Drug-related", "Weapons",
                       "Robbery", 
                       "Prior Offense", "Previous Turnout", "Non-Eligibility",
                       "Income Below", "Income Above", "Income Unk")
    
    tab.b <- data.table(t(tab[, c(1, 3)]))
    tab.b$type <- c("mean", "w_mean")
    tab.b
  }
}

results.c <- do.call('rbind', results)
results.c1 <- results.c[type == "mean"]
results.c2 <- results.c[type == "w_mean"]
results.c1$type <- results.c2$type <- NULL

round(t(tab.o), 3)

cbind(round(apply(results.c1, 2, sd), 3), round(apply(results.c2, 2, sd), 3))
